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ABSTRACT 



Aims. The question is addressed whether stellar differentially rotating radiative zones (like the solar tachocline) excite nonaxisymmet- 
ric r-modes which can be observed. To this end the hydrodynamical stability of latitudinal differential rotation is studied. The amount 
of rotational shear required for the instability is estimated in dependence of the character of radial stratification and the flow patterns 
excited by the instability are found. 

Methods. The eigenvalue equations for the nonaxisymmetric disturbances are formulated in 3D and then solved numerically. Radial 
displacements and entropy disturbances are included. The equations contain the 2D approximation of strictly horizontal displacements 
as a special limit. 

Results. The critical magnitude of the latitudinal differential rotation for onset of the instability is considerably reduced in the 3D 
theory compared to the 2D approximation. The instability requires a subadiabatic stratification. It does not exist in the bulk of con- 
vection zone with almost adiabatic stratification but may switch on near its base in the region of penetrative convection. Growth rates 
and symmetry types of the modes are computed in dependence on the rotation law parameters. The S 1 mode with its transequatorial 
toroidal vortices is predicted as the dominating instability mode. The vortices show longitudinal drift rates retrograde to the basic 
rotation which are close to that of the observed weak r-mode signatures at the solar surface. 
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1. Introduction 

The differential stellar rotation may excite other types of motion 
via instability. The possible transmission of rotational energy to 
other types of motion may be relevant for various astrophysical 
processes. Excitation of r-modes of global oscillations in differ- 
entially rotating neutron stars is considered as a source of de- 
tectable gravitational waves (Watts et al. 120031 1. Knaack et al. 
( I2005I I interpreted the large-scale structures in magnetic fields 
of the Sun as signatures of r-modes which may in turn result 
from an instability. 

The stability problem of differential rotation is also relevant 
for the dynamics of the solar tachocline (Oilman 12005 j ). The 
tachocline is the thin shell beneath the convection zone where 
the rotation pattern changes strongly. Beneath the tachocline the 
solar rotation is practically uniform. Above the tachocline, the 
rotation rate varies with latitude as observed at the solar surface. 
Inside the tachocline, a transition from differential to uniform ro- 
tation occurs with increasing depth. The question is whether this 
tachocline is hydrodynamically stable or not. If it is, the idea 
that the site of the solar dynamo is beneath the convection zone 
is hard to accept. 

The tachocline thickness is about 4% of the solar radius. The 
tachocline is located mainly if not totally beneath the base of the 
convection zone at Rin = 0.71 3/?© (Christensen-Dalsgaard et al. 
1991; Basu & Antia 1997) in the uppermost radiation zone. 

The stability/instability of the solar tachocline is also related 
closely to the lithium problem. The lithium at the surfaces of 
cools MS stars slowly decays with a characteristic time ~ 1 Gyr. 
The primordial lithium is destroyed at temperatures greater than 
2.7 Mio K which is already exceeded 42.000 km beneath the bot- 



tom of the convection zone. Evidently, the tachocline should not 
be too unstable otherwise the downward transport of the lithium 
may easily be to strong. Nevertheless, the diffusion coefficient 
for the lithium must exceed the molecular diffusion by one or 
two orders of magnitude. 

The hydrodynamical stability problem has extensively been 
studied in 2D approximation of purely toroidal disturbances. 
Symmetry types and growth rates of the 2D unstable modes 
are known (Watson 119811 Dziembowski & Kosovichev 119871 
Charbonneau et al. 1 19991 1, and also the weakly nonlinear evolu- 
tion of the instability has been described (Garaud 2001'; Cally 
2001). The 2D approximation neglects the radial displacements 
which can be expected as small in stably stratified radiative 
shells where the buoyancy frequency is much higher than the 
rotation rate, » Q (Watson [19811 1. This condition is not ful- 
filled in stellar convection zones. 

The present paper overcomes the 2D approximation by al- 
lowance for radial displacements. Poloidal motions and entropy 
disturbances are thus included. Our formulation contains the 
2D approximation as a special limit of large parameter A = 
{N/(rQ.) » 1 ; where { is the radial scale of the disturbances and 
r is the radius. We shall see that the most unstable modes have so 
small radial scales that Ail and this condition is by far not ful- 
filled. The minimum amount of differential rotation for onset of 
the instability is considerably reduced compared to the 2D case. 
More important is that the instability does not exist in the limit 
of A — > 0, so that differential rotation is stable in convection 
zones with their almost adiabatic stratification. The instability 
may, however, switch on in the region of penetrative convection 
near the base of the convection zone. Such a near-base instabil- 
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ity may be the reason for the difference in latitudinal profiles 
of angular velocity between the top and the bottom of the solar 
convection zone (Charbonneau et al. 1999). If this is the case, 
transequatorial vortices (unstable S 1 modes) should be present 
near the base. The rates of (retrograde) drift of the vortices are 
similar to that of the r-modes signatures inferred by Knaack et 
al. (120051 1 from solar magnetograms. 

2. The model 

The latitudinal dependence of the angular velocity Q. on the 
sun can be approximated by an expression including cos^ 9 and 
cos'^ 9 terms so that 

n = Oo(l -a[(l-f)cos^e + fcos*6)). (1) 

Here Q.o is the equatorial angular velocity, a is the normalized 
equator-pole difference in the rotation rate, and / is the frac- 
tion of the cos'* 9 term contribution to that differential rotation. 
The rotation is assumed slow enough not to deform the spherical 
symmetry of the star 

At the solar surface it is a ^ 0.3 and / ^ 0.5 (Howard et al. 
I1983I I. The latitudinal shear varies only slightly with depth in the 
bulk of the convection zone but it shows a characteristic change 
near its base (Charbonneau et al. 119991 ). The amplitude a{l - /) 
of the cos^ 9 term remains almost constant up to the base and 
starts decreasing in the deeper tachocline only while the fraction 
/ of cos^ contribution drops to practically zero near the base 
(cf. Fig. 10 of Charbonneau et al. il999D . 

The stratification is characterized by the buoyancy frequency 

N^-jr^, (2) 
Cp dr 

where g is the gravity, Cp is the specific heat at constant pressure, 
and s is the specific entropy. 

We address the linear stability problem with the small dis- 
turbances depending on time as exp(-ia;f). A positive imaginary 
part of the eigenvalue u) means an instability. The radial scales 
of the disturbances are assumed small compared to the stellar 
radius while the equations are global in both the horizontal di- 
mensions. The dependencies on radius and longitude are taken 
as Fourier modes exp(iOT0 + \kr). 

2.1. Equations 

The linear equations for small perturbations in differentially ro- 
tating fluids with toroidal magnetic fields are akeady given by 
Kitchatinov & Riidiger (120081 ). Here, the nonmagnetic version 
of the equations is considered in detail. The equations are for- 
mulated for normalized parameters (the rules of conversion to 
physical variables are given below). The equation for the poten- 
tial W of the toroidal flow reads 

{&) - mQ) [LW] = -i^l [LW] - m — ^—r-. W 



dj-L 



dV 



(3) 



where a> - w/Qq is the normalized eigenvalue, Q - Q/Qq is 
the normalized rotation rate, /j. - cos 9, V is the poloidal flow 
potential, 

L^^{l-^']^-^^ (4) 



is the angular part of the Laplacian operator, and 



A = 



(5) 



is the key parameter for the influence of the stratification. The 
diffusion terms are characterized by the parameters 



ey = 



(6) 



where v and;^' ^6 the microscopic viscosity and the thermal con- 
ductivity. 

Apart from the Eq. (O for the toroidal flow, the complete 
system of three equations includes the equation for poloidal flow, 

icb - mCl) [lv] = -i^ {lv) - [ls] + 

+ 2m[^^V + (l-f/)——\- 

. X , d(uCl)dW .dCl 

- l^iQ {lw) - 2(1 - - 2m^— W, (7) 



and the equation for the normaUzed entropy S , 

(oj - mCl) S = -i^S + LV. 

A^ 



(8) 



The equations (O, (|7|) and ([8]) form an eigenvalue problem 
which we solved numerically. The diffusion parameters ^ for 
the upper radiative core of the Sun are - 10"'* and ey - 2 x 
10""^ which are used in the calculations. 

The disturbances in physical units follow from their normal- 
ized values by 



iCpAT^ 
gkr 



S, P„ = [nor^/k] V, r„ = (ilor^) W. 



(9) 



The velocity field can be restored from the potentials of poloidal 
(P„) and toroidal (T,,) flows. 



er f„ eg I I dT„ d^P,, 

u - -^LPu — + —-— 

r'- r \sm0 o0 drd9 



r \ do sin 9 drd(f> 



(10) 



(Chandrasekhar 1 96 1), where e,-, eg, and are unit vectors in 
the radial, meridional and longitudinal directions. 

Without rotation (Q — > 0) and for small diffusion the equa- 
tions ([3]), (|7| and ([8]) reproduce the spectrum 



2 l{l+l)N^ 



(11) 



1 



of ^-modes. More related to the stability problem is the limit 
of very large /i-parameter ^ which leads to the following 2D 
approximation. 

2.2. 2D approximation 

The ratio of N^/Q.^ in stars can be so large (~ 10^ in the upper 
radiative core of the Sun) that A- (|5]) can also be large in spite of 
short-wave approximation in radius, kr '» 1 . In the limit of large 
A^ the above equation system reduces to its 2D approximation. 
In leading order of this parameter Eq. ^ gives S - 0. Then it 
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follows from (|8j that V = and Eq. ([3) reduces to the standard 
equation of 2D theory of Watson (119811 ). 



(w - mQj [LWj = -m — ^-^^-^ W, 



(12) 



describing toroidal flows on spherical surfaces. 

The 2D approximation is justified for stable oscillations with 
not too short radial scales so that A remains large. Its validity for 
stability problem is less certain because the radial scales of most 
rapidly growing modes are not know in advance and the value of 
kr for those modes is normally so large that Ail (Kitchatinov 
& Rudiger 2008j. 

For rigid rotation Eq. ( fT2] i provides the eigenvalue spectrum 



/ ^ 
a> — niQ. 1 

I 1(1 + 1) 



/ = 1,2,... 



(13) 



of the r-modes (Papaloizou & Pringle 119781 1. Instabilities can 
emerge with nonuniform rotation. The necessary condition for 
instability is that the second derivative, d^((l - //^)Q)/dyU^, 
changes its sign (Watson [198T). For the angular velocity profile 



2« 



the condition demands 



1 



a > 



4n + 1 



(14) 



(15) 



It should be a > 0.2 for n - I. For n - 2, i.e. the cos"* 0- 
profile ( fT4l i. the amplitude a of differential rotation for the onset 
of instability may reduce considerably, a > 1/9. This is why the 
cos** 0-term is kept in the differential rotation profile ([T). 

The profile ([T) for the Sun is, of course, an approximation. 
Higher order terms in cos^ may also be present. Reduction of 
the instability threshold due to the higher order terms is, how- 
ever, less significant and they are relevant only to the near-polar 
regions. Here, the results are presented for the rotation law ([l)- 

2.3. Symmetry types 

The eigenmodes provided by both the 2D approximation (IT2l i 
and the full 3D equation system of section ITT] possess definite 
equatorial symmetries. We use the notation Sm for the modes 
with symmetric relative to the equator potential W of toroidal 
flow and the notation Am for antisymmetric W (m is the az- 
imuthal wave number). The symmetry convention is same as 
used before (Charbonneau et al. 1999; Kitchatinov & Riidiger 
[2008). Note that the eigenmodes combine W of definite symme- 
try with S and V of opposite symmetry, i.e., S and V are sym- 
metric for Am modes and antisymmetric for Sm. The velocity 
field for Am modes has antisymmetric ug and symmetric Ur and 
about the equatorial plane, and the other way round for Sm 
modes. 

Watson ( 119811 ) proved that only nonaxisymmetric modes 
with m = 1 and m - 2 can be unstable in the 2D approxima- 
tion. Our 3D computations yield the same conclusion. 

3. Results 

Figure [T] shows the critical shear amplitudes a as functions of A 
for / = 0. This is the case considered by Watson. As it must be, 
his results are reproduced in the limit of large A. Note, however, 
that the most easily excited modes have ^ < 1 . 
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Fig. 1. Neutral stability lines for / = in rotation law ([T]l. The 
instability region is above the lines. Only Al and S2 modes are 
unstable. The lines approach the marginal a-values of Watson 
theory for large A. Note that the most unstable modes have A < I. 
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Fig. 2. The same as in Fig.[T]but for / = 0.5. The critical shear 
for onset of the instability is reduced, the newly appeared unsta- 
ble modes S 1 and A2 are excited most easily. 



Even for large NjQ., the small radial displacements are sig- 
nificant for the instability. The reason is that the most unstable 
modes have short radial scales. It can be seen from ( fTOl i that zero 
radial velocity would exclude the whole class of poloidal distur- 
bances. The ratio of horizontal (mh) to radial velocities in a cell 
of poloidal flow of different radial {() and horizontal (H) scales 
can be estimated as u^lu^ ~ Hfl. Horizontal velocity of poloidal 
flow can thus remain important in spite of small Ur if the ra- 
dial scale is much shorter than the horizontal one. The poloidal 
(interchange-type) disturbances are so significant for the insta- 
bility that the critical latitudinal shear for onset of the instability 
reduces from a - 0.28 (value of Watson) to a = 0.21. However, 
this smaller value still ensures the tachocline to be stable. 
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Fig. 3. Isolines for the normalized growth rates cr - 3(w)/f2o of 
most rapidly growing modes. Dominating symmetry types are 
indicated. Dashed lines show the watershed between different 
symmetry types. 



We find that the instability disappears when the stratification 
approaches adiabaticity (A^ 0). The instability of differential 
rotation does thus not exist in convection zones. 

Fig. |2] re veals a striking destabilizing effect of finite / in the 
rotation law ([T]). The threshold value of a for marginal stability 
is further reduced compared to the case of / = 0. New unstable 
modes appear and the S 1 mode is now preferred. 

We find the instability as rather sensitive to the details of the 
rotation law. The growth rates of the unstable modes in depen- 
dence on a and / are shown in Fig. |3] The length scale A was 
varied to find the maximum growth rates shown in the plot. The 
dashed line separates the regions of different symmetry types. 
Surprisingly, even the symmetries of the most rapidly growing 
modes depend on the shape of the rotation law. The shape of the 
rotation law is the result of the interaction of the turbulence and 
the basic rotation in the solar/stellar convection zone. 

3.t. Angular momentum transport 

The depth dependence of the solar rotation law known from he- 
lioseismology can be interpreted in light of the presented results. 
The rotation law in the bulk of the almost adiabatically stratified 
convection zone is stable. In the region of penetrative convection 
near the base the stratification changes to subadiabaticity and the 
instability can exist. If it exists, then it reacts back on the differ- 
ential rotation to change it towards a stable profile with f - 0- 
Our linear computations cannot describe this nonlinear process 
but they can probe for the sense of angular momentum trans- 
port. Figure |4] shows that the instability indeed tends to reduce 
the differential rotation. The plot shows the angular momentum 
flux {ugu^} after longitude-averaging as a function of the latitude. 
The coiTelation is negative (positive) for the northern (southern) 
hemisphere. The angular momentum is thus transported from the 
equator to the poles. The plot was constructed for slightly super- 
critical S 1 mode which should be active if the differential rota- 
tion is reduced to the marginally stable value (Fig.O. 




Fig. 4. Meridional flux of angular momentum for slightly super- 
critical (cr - 0.001) SI -mode for / - 0.2 with normalized ve- 
locities. 
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Fig. 5. Streamlines of toroidal flow for the same mode as in 
Fig- El Full and dotted lines show opposite senses of circulation. 



3.2. The flow pattern 

The streamlines of the toroidal flow for the same S 1 mode are 
shown in Fig.|5] Close to the equator the flow represents transe- 
quatorial vortices. The flow pattern drifts in longitude against 
the direction of rotation in the corotating frame with rates shown 
in Fig. |6] If the deep-seated vortices were observable (e.g., due 
to disturbance of the large-scale magnetic field), the observer 
would see the frequencies 



V = m (vq - ve) + Vq V, 



(16) 



where Vq - 450 nHz is the equatorial rotation frequency of the 
base of the convection zone, ve = 3 1 .7 nHz is the orbital rotation 
frequency of the Earth, and v is normalized drift rate of Fig. |6] 
The equation gives the range of frequencies from about 120 to 
160 nHz for the range of drift rates of Fig. |6] The frequencies 
correspond to periods from about 120 to 160 days for the method 
of analysis of synoptic maps of solar magnetic fields by Knaack 
et al. ( I2005I I. Two of the periods which Knaack et al. interprete 
as signatures of r-modes fall into this range. 

4. Discussion 

Rotation laws in stably stratified stellar interiors with strong 
enough latitudinal gradient are hydrodynamically unstable 
against nonaxisymmetric disturbances or - what is the same - 
the r-modes are excited for large enough differential rotation. As 
their drift rate is retrograde with an amplitude of about 10% of 
the rotation rate they should be observable. 
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Fig. 6. Normalized drift rates %(cj)/Q.o - 1 in corotating frame 
for slightly supercritical (cr - 0.001) Sl-modes as function of 
/-parameter of rotation law ([T]|. 
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They are only excited, however, in subadiabatically stratified 
radiative zones and convective overshoot regions. Hence, their 
existence directly indicates the extended regions of latitudinal 
shear in the deep stellar interior beneath the proper convection 
zone. If this region - like the solar tachocline - is thin the in- 
tensity of the r-mode oscillation should be much smaller than 
in stars without tachoclines. New calculations for rotation laws 
with radial profiles of the angular velocity are necessary to de- 
velop this new r-mode seismology. In the present paper the radial 
shear of rotation is not included that can be justified only if the 
angular velocity varies in radius on a scale larger than the radial 
wavelength £lr/N of unstable excitations. 

The critical latitudinal shear for the excitation of the modes is 
not very small. The simplest theory without radial perturbations 
and with a simplified parabolic rotation law yields a critical lat- 
itudinal shear of 28%. We have shown with an improved math- 
ematics that the real value is less. It is reduced to 21% for the 
same rotation law but with a 3D theory. The critical shear value 
is further reduced if the rotation law contains a higher order 
term of cos'* 0. Nevertheless, the critical shear rate remains larger 
than (say) 10%. A rotation law with a slightly smaller latitudi- 
nal shear (driven by the turbulence in the convection zone) could 
stably exist in the stellar interior without any decay. We know, 
however, that the solar core rotates almost rigidly. This is only 
understandable with the idea that the slender solar tachocline is 
due to an extra effect, e.g. by the Maxwell stress of large-scale 
magnetic fields. They may be of fossil origin as their amplitudes 
need not exceed (say) 1 Gauss (cf. Rudiger & Kitchatinov |2007l l. 

If this magnetic concept for the solar tachocline is true and 
if the tachocline is indeed hydrodynamically stable then also the 
very slow decay of the observed lithium abundance becomes un- 
derstandable. The slight increase of the diffusion coefficient by 
one or two orders of magnitude can easily be explained with 
slow horizontal motions of order cm/s (see Riidiger & Pipin 
1200 11 1 or by radial plumes overshooting from the convection 
zone (Blocker et al. |19981 l. 
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